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The density functional tight binding approach (DFTB) is well adapted for the study of point 
and line defects in graphene based systems. After briefly reviewing the use of DFTB in this area, 
we present a comparative study of defect structures, energies and dynamics between DFTB results 
obtained using the dftb+ code, and density functional results using the localised Gaussian orbital 
code, AIMPRO. DFTB accurately reproduces structures and energies for a range of point defect 
structures such as vacancies and Stone- Wales defects in graphene, as well as various unfunctionalised 
and hydroxylated graphene sheet edges. Migration barriers for the vacancy and Stone-Wales defect 
formation barriers are accurately reproduced using a nudged elastic band approach. Finally we 
explore the potential for dynamic defect simulations using DFTB, taking as an example electron 
irradiation damage in graphene. 



I. INTRODUCTION 

The family of carbon structures is much larger than its 
most notable components and the number of new mem- 
bers synthesized each year makes it hard to categorize all 
carbon forms. ^ Our ability to describe computationally 
the structure of realistic carbon systems faces the ad- 
ditional difficulty represented by the presence of native 
defects, which often dominate the mechanical, electronic 
and chemical properties of their host material. Further- 
more, defect combinations can also serve as elemental 
topological transformations that when applied to origi- 
nal perfect forms generate more complex structures. 

Modern electron microscopes can access subnanomet- 
ric spatial resolutions and thus are now able to image in- 
dividual defects in nanostructures. However image inter- 
pretation is not straightforward and accurate structural 
models are still required for an in-depth understanding 
of defective atomic structures. In this context computa- 
tional simulations represent a necessary complementary 
tool both for the interpretation of experimental data and 
for a deeper understanding of the specific physical and 
chemical properties at defective sites. 

Classical density functional theory (DFT) methods 
have been shown to describe with a high accuracy the 
structure and electronics of defects in carbon materials. 
Point defects can be simulated using either clusters or 
periodic structures formed from hundred of atoms. How- 
ever more extended defects such as extended dislocation 
lines^, turbostratic or misaligned graphite"^ or amorphous 
carbon^ require significantly larger models. Until re- 
cently density functional techniques have been limited 
to carbon models containing only few hundreds of atoms. 
This limitation is going to be overcome by very recent im- 
provements in filtration techniques on Gaussian basis sets 
that can massively reduce computational time and mem- 
ory requirements.^ Thanks to these methodological de- 
velopments combined with the continuous speeding up of 
computing systems, routine study by full DFT of atomic 



structures with many thousands of atoms is becoming a 
realistic task. 

However, beside ground state determination, the study 
of the dynamics of defective crystals remains too com- 
putationally intensive to tackle within the framework of 
current density functional theory. Examples of these spe- 
cific problems include the growth mechanism of carbon 
nanostructures, their thermal or mechanical stability, the 
effects of structural reorganization induced by high en- 
ergy particle irradiation, and the diffusion of extended de- 
fects such as vacancy clusters or dislocations. This range 
of problems are generally investigated using molecular 
dynamics or Monte Carlo techniques where forces and 
energies are evaluated using empirical or semi-empirical 
approachesi^ However many dynamical problems in car- 
bon nanostructures involve mechanisms of bond break- 
ing and reconstruction that are usually poorly estimated 
by computational methods parametrized on ground state 
configurations. 

The usage of density functional tight binding (DFTB) 
for the study of dynamics and reorganization of com- 
plex carbon structures represents an extremely helpful 
compromise between accuracy and speed. ^^^'^ Density 
functional tight binding parameters, in particular those 
derived for carbon,^! are highly transferable, overcom- 
ing the main limitations of empirical and semi-empirical 
techniques at the reduced computational cost of stan- 
dard tight binding. This is demonstrated by a wide range 
of problematics covered by a number of studies on car- 
bon based systems. DFTB has been employed for in- 
stance in the study of the structure and energetics of 
point defects in single walled carbon nanotubesii"— and 
more extended defective structures such as screw dislo- 
cations in multi walled carbon nanotubes^^ and bond- 
ing between fuUerenes and nanocones.^^ The capability 
of carbon DFTB parameters to reproduce complex reby- 
dridization phenomena has been shown for example in 
the simulation of the thermal induced graphitization of 
nano-diamond surfacesJ^ in monoatomic carbon chains 
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formation at axial strain applied to carbon nanotubes,^9 
and carbon nanostructure growth through C2 additionf2i 
DFTB molecular dynamics has been also employed in the 
study of self-assembling processes of fullerene cagea^Sri^i 
and nanotubes.^^ 

In order to justify the use of DFTB approaches it is 
necessary to benchmark the results against more conven- 
tional approaches such as full DFT calculations. Surpris- 
ingly in the literature there is a lack of detailed compara- 
tive studies of this nature for defective carbon nanomate- 
rials and for this reason we have undertaken the current 
work. We have first chosen a range of well characterised 
intrinsic point defect structures in graphene in order to 
benchmark the DFTB optimised geometries and forma- 
tion energies. This is followed by a study of more com- 
plex unterminated and functionalised edges to quantify 
its capability with extended structural defects. As well 
as static ground state structures it is important to anal- 
yse the performance of the DFTB approach further from 
equilibrium and for this reason we next study defect diffu- 
sion and formation/annihilation barriers. Finally we take 
as a complex dynamic example the formation of points 
defects in carbon nanomaterials under the influence of 
an electron beam. As we show, DFTB is capable of re- 
markably accurate reproduction of full DFT calculations 
at a fraction of the computational cost, justifying its use 
in a wide range of structurally complex carbon nanoscale 
problems. 



II. COMPUTATIONAL METHOD 

DFT calculations are conducted using the AIMPRO2.0 
code2&i21 and the DFTB approximation as implemented 
in the dftb-|- code^^, using comparable cells and k-point 
meshes. The AIMPRO calculations are performed under 
the local density approximation using a localised Gaus- 
sian basis set with 22 independent functions per car- 
bon atom (12 per hydrogen and 40 per oxygen). Finite 
temperature Fermi smearing is used to control electron 
state population near to the Fermi level with tempera- 
ture kT=0.04 eV. Spin polarised calculation have been 
performed for open shell configurations (carbon mono- 
vacancy, zigzag and Klein edges). Pseudopotentials are 
taken from Hartwingser-Goedecker-Hutteri^ DFTB pa- 
rameters have been derived by M. Elstner et al.^° Both 
the DFT and DFTB calculations are fully self-consistent. 
No additional functionality such as Van der Waals cor- 
rections are used within the DFTB calculations. 



III. POINT AND LINE DEFECTS IN CARBON 
NANOSTRUCTURES 

In this section we present a comparative study on the 
structure and formation energies of topological defects 
in graphene and grapheme structures. We take as our 
first test system a series of standard intrinsic point de- 



fect structures in the graphene lattice, namely a single 
vacancy, a 5-8-5 divacancy pair, a "Stone- Wales" defect 
(a rotation of two carbon atoms through 90° about their 
bond centre), and an "inverse Stone- Wales" defect (ad- 
dition of a C2 pair to the graphene lattice across a single 
hexagon). Between them these defects contain a vari- 
ety of local bonding including undercoordinated carbon 
atoms, dilated bonds, local out-of-plane distortions and 
combinations of resonant and localised single-/double- 
bonds. As such they represent a stringent test for DFTB. 
Through the use of an infinite graphene nanoribbon we 
then examine extended defects in the form of untermi- 
nated graphene sheet edges. These once again exhibit a 
range of bonding states including localised triple bond 
character (armchair and reconstructed zigzag edges), ex- 
tended metallic states (zigzag edges) and singly coordi- 
nated carbon atoms (Klein edge). 

In figure [T] we present structural models for the dif- 
ferent types of point defects in graphene and several 
graphene edge configurations. In the figure we report the 
most notable bond lengths as obtained after optimization 
with DFT (blue values) and DFTB in its self consistent 
charge formulation (red values). The DFT results are 
detailed further in Rei^. 

The DFTB structures are in excellent agreement with 
the DFT data for all carbon coordinations, from the 
single coordination of a carbon atom at a Klein edge, 
through the double coordinate carbons at the zigzag and 
armchair edge or close to the vacancy site, to the modi- 
fied triple coordinated states as in the Stone- Wales and 
inverse Stone- Wales defect (C2 addition). All DFTB- 
DFT bond lenght discrepancies are lower then 4% and 
most around 2%. 

Formation energies are presented in table U and once 
again agreement is excellent, with DFTB point defect 
formation energies matching DFT values to within 1.5% 
(and most less than 1%). The largest error is in the Klein 
edge formation energy (19.7%), which is understandable 
given that this is a physically unstable edge structure 
(repetition of the unit cell and breakage of the symmetry 
leads to spontaneous pairwise rebonding of the underco- 
ordinated atomsSS). The other edge formation energies 
deviate from the DFT values by 9.7%, 1.8% and 1.8% 
respectively. 

Pristine graphene edges have dangling bonds at the 
edge atoms. Rehybridization (as considered above) or 
H-termination is the simplest way to saturate these dan- 
gling bonds, while inducing only small edge strain. Edge 
functionalization by halogens or more complex functional 
groups (-0H or -SH) has been shown to induce a sig- 
nificant strain along the ribbon edge, through steric 
hindrance, electrostatic repulsion between groups, inter- 
group bonding, etc. Being energetically unfavourable, 
this strain can be relieved via out-of-plane distor- 
tions. Specifically, hydroxyl (-0H) terminated graphene 
nanoribbons of different widths have been shown to com- 
pensate the induced strain by forming a localised out-of- 
plane static ripple along the graphene sheet edge.'^*' A key 
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Figure 1: Structure of intrinsic point defects and unterminated edges in graphene. Values in the figures represent bond lengths 
in Angstroms obtained from DFT (blue) and DFTB (red). DFT values from Reference.— 



Point defects (eV) 


DFT DFTB 


Mono- vacancy 


7.40 7.51 


5-8-5 Di vacancy 


8.25 8.19 


Stone- Wales 


4.86 4.85 


Inverse Stone- Wales 


6.37 6.40 


Edges (eV/A) 


DFT DFTB 


Zigzag edge 


1.34 1.21 


Armchair edge 


1.10 1.08 


Klein edge 


2.22 1.78 



Zigzag reconstructed edge 1.09 1.07 

-OH terminated armchair ribbon -2.26 -3.53 



Table I: Formation energies for point defects and edges in 
graphene obtained by DFT and DFTB. DFT edge values from 
References2^i^. 



consequence of these functionalised nanoribbon edges is 
that both electronic and mechanical properties can be 
tuned,^ 

In figure [H we present the structural model of an hy- 
droxyl terminated graphene nanoribbon. Both DFT and 
DFTB calculations confirm that a rippled configuration 
is more stable than any flat structure, allowing the ribbon 
edge to relieve strain through out of plane distortioni^ 

Once again, bond lengths obtained with DFT and 
DFTB generally correlate well. The only discrepancy 
occurs at hydrogen bridges whose lengths are over- 
estimated by DFTB. This bond dilation is also refiected 
in the edge formation energy in table U where the DFTB 
edge formation energy is too energetically stable com- 
pared to the DFT result. DFTB is known to tend towards 



overbinding for hydrogen-X bonds, and the dftb-|- code 
includes a damping correction to the short range contri- 
bution to the sec interaction for hydroge n^^i'^^ which 
we did not use here. In addition we can also not ex- 
clude the effect of the limited size of the DFTB basis 
or the possibility of an incorrect estimation of the -OH 
group chemical potential, and further studies are needed 
to fully explain this difference. 

In general we find a good correspondence in the ener- 
getics (formation energies) and structural characteristics 
obtained by DFTB with DFT resuhs, for both intrin- 
sic point and line defects in graphene, with reasonable 
structural and energetic agreement for heteroatom sys- 
tems given the limitations discussed above. These re- 
sults support the use of DFTB in future, notably for 
problems which can not be treated at the DFT level due 
to their size and complexity, or due to the necessity for 
long trajectories, for example DFTB-MD simulations on 
dynamics of graphene ribbon rippling (propagating and 
stationary waves along the edge). 



IV. DEFECT DYNAMICS: FORMATION AND 
MIGRATION BARRIERS 

The description of complex dynamical phenomena such 
as the reorganization of carbon material when exposed 
to thermal treatments can be decomposed into several 
elementary processes such as atom and vacancy migra- 
tion and nucleation, bond rotation, and carbon atom 
re-hybridization. Therefore, it is fundamental to derive 
accurately minimum energy reaction paths and activa- 
tion energies of these elemental transformations in order 
to describe precisely more global transformations. How- 
ever, the ability of carbon to re-hybridize can render even 
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Figure 2: Structure of -OH terminated armcliair graphene 
nanoribbon where values in the figures represent bond length 
(in A) obtained from DFT (blue) and DFTB (red). DFT 
values from Reference^. 



Barrier (eV) 




DFT 


DFTB 


Vacancy Migration 


1.37^ 


1.29 


Stone- Wales Formation 




10.4 


Stone- Wales Annihilation 


4.4^* 


4.7 



Table II: Calculated barriers (eV) for point defect formation 
and migration in graphene obtained by DFT and DFTB using 
the nudged elastic band method. 

elemental reaction paths highly complicated, often intro- 
ducing intermediate metastable configurations. 

The nudged elastic band (NEB) method^^i^ is a useful 
techniques for exploring in an efficient and automatic way 
a large region of the configuration space and derive com- 
plex minimum energy paths. However, a large number 
of intermediate images is generally required for obtain- 
ing with a good accuracy the saddle point configuration 
and associated activation energies. Furthemore the NEB 
method converges usually only after a large number of 
optimization steps. The high computational cost of the 
NEB technique represents a strong limit for its usage in 
the framework of high level computational approaches. 

A good compromise can be obtained using the den- 
sity functional tight binding theory (DFTB). A combined 
NEB-DFTB approach has already been shown to give re- 
sults comparable to other higher-end techniques'^. In the 
context of layered materials, DFTB-NEB has been suc- 
cessfully employed in the study of single vacancies and va- 
cancy complex migration in boron nitride monolayers f2S. 
Here we compare the possibilities of a DFTB-NEB ap- 
proach with equivalent more time consuming DFT-NEB 
calculations for the study of topological transformations 
in carbon. Two example are presented where bond break- 
ing occurs: mono-vacancy migration and Stone- Wales 
bond rotation in a graphene plane. 

Mono-vacancy diffusion in graphene occurs when a 
doubly coordinated carbon close to a vacancy site breaks 
its two covalent bonds for rebonding with the opposite 



atom pair neighboring the vacancy (the complete migra- 
tion path also involves a bonding rearrangement around 
the vacancy core but this has a very low activation bar- 
rier). The minimum energy path saddle point obtained 
by DFTB-NEB corresponds to a configuration where the 
migrating carbon atom lies in the middle point between 
its initial and final position. The graphene sheet under- 
goes a slight asymmetric out-of-plane deformation allow- 
ing the migrating atom to locate itself at the center of 
a compressed tetrahedron. The DFTB activation energy 
we obtain is 1.29 eV, which is in an extremely good agree- 
ment with the value of 1.37 eV obtained by an analogous 
DFT-NEB study37 (6% underestimation). 

The 90 degree bond rotation required for the forma- 
tion of a Stone- Wales defect involves the breaking and 
reconstruction of two covalent bonds. Using DFTB-NEB 
we obtain an activation energy for defect formation of 
10.4 eV, with the corresponding annihilation barrier for 
the reverse reaction of 4.7 eV. These compare extremely 
well to equivalent DFT values of 9.2 eV and 4.4 eV 
respectively^-, with DFTB thus overestimating the DFT 
values by 13 and 7% respectively. 

These activation barrier calculations represent a strin- 
gent test for DFTB, passing through structures which are 
far from the equilibrium defect ground states. The ex- 
cellent agreement between DFT and DFTB justifies the 
use of DFTB in calculations of dynamic systems, and in 
the following section we show an example of this where 
we use DFTB to examine atom loss during electron irra- 
diation of carbon nanostructures. 



V. DEFECT DYNAMICS: ELECTRON 
IRRADIATION IN CARBON 
NANOSTRUCTURES 

Electron irradiation is an unavoidable and generally 
unwanted side effect when high energy electrons are used 
for imaging and analysis of nanostructures, such as in 
transmission electron microscopy (TEM), but it can also 
be used to deliberately restructure a carbon nanosystem 
in order to tune its mechanical and electronic properties. 
In this context it is desirable to be able to finely describe 
the probability that a specific structural transformation 
occurs under electron irradiation and how this probabil- 
ity depends on the energy of the electron beam. 

Electron irradiation effects in carbon materials can be 
explained mostly by direct elastic scattering between the 
relativistic electrons of the beam and the atomic nucleus. 
For a given electron beam energy atoms can only be sput- 
tered along directions for which the transfered energy 
is above a certain emission direction dependent energy 
threshold. 

An analytic expression for the differential cross sec- 
tion as a function of the emission direction has been 
derived by Motli^ and a useful approximation of the 
original expression has been obtained by McKinsley and 
Feshbach42. Total emission cross sections are derived by 
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integrating the differential cross section over the sohd an- 
gle defined by the possible emission directions at a given 
electron beam energy. 

In Ref4i we proposed a methodology for deriving 
anisotropic emission energy threshold maps using ex- 
tended molecular dynamics simulations. The procedure 
consists in imparting an initial momentum to the atom 
to sputter (direction and speed). The system is succes- 
sively allowed to evolve in a microcanonical ensemble and 
at the end of the simulation the final position of the atom 
is analyzed. The MD simulation is repeated increasing 
progressively the initial speed, up to the critical limit for 
which the atom is ejected. The procedure is reiterated 
for a number of different emission directions. 

To obtain accurate energy threshold maps a high num- 
ber of directions have to be considered and the step size 
used to increase the initial velocity should be sufficiently 
small. In the case of perfect graphene a map is ob- 
tained performing about 10000 molecular dynamics cal- 
culations, each of 200 MD steps on a structure containing 
200 atoms. This number of calculations is too computa- 
tionally demanding for standard density functional tech- 
niques, but DFTB can produce a full emission map at an 
affordable computational cost. This techniques has thus 
be successfully used in the study of sputtering in perfect 
and defective graphene and monolayer boron-nitride as 
well as in the study of irradiation induced bond rotations 
in graphene^""—. 

In figure [3] we present the DFTB-MD derived emis- 
sion energy threshold map for a carbon atom from a 
graphene plaine. DFTB estimates the minimum ejection 
energy to be around 23 eV, corresponding to an emission 
direction orthogonal to the plane. A recent work has 
compared the DFT and DFTB emission energies obtain- 
ing, in the DFT case, a value of 22.2 eVM The excellent 
agreement between the DFT and DFTB values (less than 
4% difference) validates the use of DFTB for sputtering 
simulations in carbon materials. However in the case 
of boron and nitrogen sputtering from a BN plane the 
DFTB values are lower than the DFT. This discrepancy 
has been discussed by the author as an inadequate de- 
scription of charge transfer in DFTB calculations for the 
BN system4^ 

Considering the kinematic of the scattering problem 
carbon atoms can be sputtered from a graphene plane 
by electrons with an energy above 113 keV. Experimental 
TEM studies find an electron beam energy limit situated 
between 90 and 100 keV. This discrepancy can be rea- 
sonably attributed to the well known dissociation energy 
overestimation occurring in DFT-LDA that also affects 
the DFTB parameters. This error can be considered as 
systematic and the theoretical results be can corrected 
by recalibrating on the experimental values. 

In Fig. |3]we present the total knock-on cross section as 
a function of the electron beam energy for a carbon atom 
in a graphene plane. A commonly used approximation 
considers the emission energy threshold as independent 
of the emission direction, an analytic expression for the 




Figure 3: Emission energy threshold map for a carbon atom 
in a graphene plane. The spherical coordinate represent the 
emission direction and the color scale the minimum emission 
kinetic energy. The equatorial plane corresponds to emission 
directions on the graphene sheet, poles to directions orthogo- 
nal to the sheet. 
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Energy (keV) 

Figure 4: Sputtering cross sections for a carbon atom in a 
graphene plane using the isotropic emission threshold approx- 
imation and the anisotropic approximation when the electron 
beam direction is orthogonal or parallel to the plane. Sput- 
tering cross sections for a carbon atom neighbouring a pre- 
existing vacancy in graphene (beam direction orthogonal) 



total cross section can then be found. This assumption 
brings however, as shown in Fig. SJ to overestimate the 
integration region and thus the emission cross sections 
and it cannot take into account the strong variation of 
cross section as a function of the electron beam orienta- 
tion in respect to the graphene. In figure S] we present 
also the cross section for a carbon atom neighboring a 
pre-existing vacancy site. The reduced coordination of 
the knocked carbon atom makes sputtering more prob- 
able than for an atom in perfect graphene. This higher 
cross section explains the vacancy clustering in graphene 
observed by TEM. 

Calculated cross sections have been used to optimize ir- 
radiation conditions in a scanning transmission electron 
microscope for reshape individual single walled carbon 
nanotubes at a nanometrical scale. ''^ Cross section val- 
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ues can also be used in kinetic Monte Carlo simulations 
of the global transformation of carbon structures under 
electron irradiation: calculated cross sections attribute a 
sputtering probability to each carbon atom, DFTB can 
successively be employed to optimize the system after 
each ejection event. 

VI. CONCLUSIONS 

The versatility of carbon in its range of bonding leads 
to a rich variety of low symmetry materials, structures 
and defects. However the corresponding size of the resul- 
tant calculations, the range of minima to explore, and 
the complex energy surfaces in non-equilibrium situa- 
tions render many such problems outside the scope of 
conventional density functional approaches. We show in 
the current study that the density functional tight bind- 
ing approach is able to successfully reproduce, with high 



quantitative accuracy, both structural and energetic data 
from full density functional calculations at a fraction of 
the computational cost. Calculations of atom knock-on 
cross sections under electron irradiation provide an ex- 
ample where DFTB calculations are able to advance and 
guide our manipulation of carbon materials at the atomic 
scale. These results confirm that DFTB represents a 
powerful tool for computationally intensive studies of car- 
bon nanomaterials. 
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